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Abstract 


The vestibulo-ocular reflex (VOR) is a well-known dual mode bifurcating sys- 
tem that consists of slow and fast modes associated with nystagmus and sac- 
cade, respectively. Estimation of continuous-time parameters of nystagmus and 
saccade models arc known to be sensitive to estimation methodology, noise and 
sampling rate. The stable and accurate estimation of these parameters are criti- 
cal for accurate disease modelling, clinical diagnosis, robotic control strategies, 
mission planning for space exploration and pilot safety, etc. 

This paper presents a novel indirect system identification method for the es- 
timation of continuous-time parameters of VOR employing standardised least- 
squares with dual sampling rates in a sparse structure. This approach permits 
the stable and simultaneous estimation of both nystagmus and saccade data. 
The efficacy of this approach is demonstrated via simulation of a continuous- 
time model of VOR with typical parameters found in clinical studies and in the 
presence of output additive noise. 


1 Introduction 

The eyes play a critical role in maintaining balance. They arc directly con- 
nected to organs of equilibrium, most importantly the inner car. Paired struc- 
tures called semicircular canals behind the cars sense motion and relay infor- 
mation to balance control centres in the brain. Health of the semicircular canals 
can be monitored by studying the vestibulo-ocular reflex (VOR). The VOR is a 
reflexive eye movement that stabilises images on the retina during head move- 
ment. Ocular responses to head perturbations consist of intermingled segments 
classified as “slow” (nystagmus) or “fast” (saccade), according to their average 
speed characteristics. Changes in VOR activity can be an indication of seri- 
ous brain or head trauma, which can negatively impact an astronaut or pilot’s 
ability to successfully achieve critical mission objectives. 

NASA’s Advanced Exploration Systems and Human Research Programs 
arc critically interested in pioneering new capabilities allowing future human 
missions beyond Earth orbit [1,2], The objective quantification of mechanisms 
that control VOR will enhance the capabilities of NASA’s human exploration 
missions by developing techniques that permit prediction of spatial disorien- 
tation and development of individualised countermeasures for astronaut crews 
that must perform critical space operations under varied gravito-inertial condi- 
tions such as launch, landing and orbital maneuvering [3]. 

It has been shown that the VOR can be accurately modelled by the NAR- 
MAX (Non-linear AutoRegressive, Moving Average exogenous) structure [4], 
Moreover, a modified extended least-squares (MELS) algorithm was devel- 
oped to estimate unbiased parameter values of non-linear Hammerstein struc- 
ture multimode systems. This two-step process adds overhead to data analysis, 
which can be significant for sufficiently long data records and batch processing. 
However, it can be reduced to one step by exploiting the nature and efficiency 
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of sparse matrices. 

Recently, it was demonstrated that the utilisation of a sparse matrix ap- 
proach is ideally suited for system identification of switched systems with ap- 
plication to the VOR [5]. This study demonstrated that a sparse matrix ap- 
proach is an efficient technique for the simultaneous analysis of nystagmus and 
saccade dynamics. However, several questions remain unsolved in this mod- 
elling and identification approach. 

The problem of continuous-time system identification from sampled input- 
output data can be divided into two broad approaches: (i) indirect methods, 
where a discrete-time model is estimated from sampled data; then an equiv- 
alent continuous-time model is calculated and (ii) direct methods, where a 
continuous-time model is obtained directly without going through the inter- 
mediate step of first determining a discrete-time model; based on concepts 
of approximate numerical integration to recreate time-derivatives needed in 
continuous-time formulations [6] . 

Here, we propose to investigate the problem of continuous-time system 
identification by developing methodology using an indirect approach. One such 
approach relies on matrix preconditioning techniques, such as standardised 
least-squares, to improve the spectral properties of the regressor matrix [7, 8], 
Often when the clustered spectrum is away from zero it results in rapid and 
robust convergence, especially when the preconditioned matrix is close to nor- 
mal. We deem this will provide more robust solutions and greater consistency 
for non-linear biological systems by circumventing issues of implementing nu- 
merical derivatives and filter selection required by direct techniques, which 
arc more challenging in a non-linear framework. Here, we focus on one crit- 
ical issue, namely, that of mapping the underlying continuous-time system to 
discrete-time for estimation, then inverse mapping back to continuous-time to 
provide physiological relevance and insight. 

The identification of continuous-time VOR par ameters is further challenged 
due to the slow and fast speed characteristics of the system. This corresponds 
to one pole close to the j co-axis in continuous-time or the unit circle in discrete- 
time whilst the other pole is located significantly inside the left hand plane/unit 
circle. To avoid aliasing of the fast mode, the system is typically oversampled 5 
to 10 times the highest known dynamics. This results in the slow mode moving 
close to the unit circle as the sampling interval approaches zero, leading to nu- 
merical instability for parameter estimation. To avoid the numerical instability 
problem we propose using a dual sampling rate technique and, thus, permitting 
stable and robust continuous-time parameter estimation. 

The organisation of this paper is as follows. In §2 we formulate the iden- 
tification problem addressed here. Section 3 introduces a dual sampling rate, 
standardised least-squares method in a sparse structure, which permits stable 
indirect estimation of continuous-time VOR parameters. This sparse matrix 
formulation is capable of simultaneously quantifying both nystagmus and sac- 
cade data. Section 4 provides results of the proposed algorithm on a simulated 
VOR model whilst §5 provides a discussion of our findings. Section 6 sum- 
marises the conclusions of our study. 
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2 Problem Statement 


Consider the continuous-time VOR system illustrated in Fig. 1 and charac- 
terised by Eqn. 1 as 



Figure 1 . System describing vestibular dynamics. Switch position 5 j: nystag- 
mus. Switch position S2: saccade 
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where F; (5) and F2 (j) are first-order dynamics for slow and fast-phase modes, 
y 1 i and y 2 ( represent the initial values at each switch time and q. r arc the num- 
ber of switches [ 4 , 5 ]. It has been demonstrated that a NARMAX description 
of VOR slow and fast-phases is [ 4 ] 
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where 8 is the Kronecker impulse function, j,k are the lags on the <5i,th and 
§2/th impulse and q. r are the number of data segments of sub-system one and 
two, respectively. The unknown model parameters arc compactly represented 
as 01,2, = [01,0! 02,02 03,03 04,04 05,05 ■ ■ ■ K h X t } T . 

Table 1 illustrates the relationship of the discrete-time (DT) parameters in 
Eqn. 2 to the underlying continuous-time (CT) parameters in Eqn. 1 and their 
inverse. Notice that the estimated overall gain is a product of the linear system 
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Table 1 . (Left) Forward and (Right) inverse relationship of NARMAX model 
parameters to underlying continuous-time parameters. 

and static non-linearity gain. Clearly, the system and its parameter mapping 
in the forward and inverse direction are non-linear and governed by z = {^27=7} 

and s = y pqqj where T is the sampling interval. 

Given this framework, the goal is to estimate the NARMAX model param- 
eters as 


min \ II (Zi, 2 - 01,20 1 , 2 ) ||; (3) 

01,2 ^ 

where Z1.2 G M. Nx 1 is a vector of outputs, 0 { 2 G R' Vxp is a matrix of regressors 
and 0 1 , 2 , G M pxl is a vector of unknown coefficients. Using the estimate 0 1,2 
and the relationship on the RHS of Table 1 we calculate the continuous-time 
VOR parameters. This approach is known as indirect system identification of 
continuous-time parameters. 

In addition, notice the expression in Eqn. 3 is a non-sparse matrix approach, 
requiring both fast and slow-phase to be analysed separately. However, using a 
sparse matrix approach by constructing a diagonal block-oriented data matrix 
and exploiting its natural sparseness as 
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it is possible to simultaneously estimate both nystagmus and saccade parame- 
ters. 
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3 Standardised Least-Squares With Dual Sampling 


In a clinical setting the VOR is appropriately sampled taking into consideration 
fast-phase dynamics. However, a suitable sampling rate for saccade results in 
slow-phase dynamics being highly oversampled. Although it is desirable to 
oversample system dynamics to avoid aliasing, too high a sampling rate can 
lead to numerical instability in system identification resulting in biased esti- 
mates [9, 10]. To avoid numerical instability we propose using a dual sampling 
rate approach to analyse both modes of VOR. This dual sampling rate method 
is easily achieved by modifying the sparse structure in Eqn. 4 as 


" z il _ i>* oi rer 

Z 2 0 02 


(5) 


where the superscript denotes that the nystagmus signal has been down 
sampled to achieve better numerical results. Here, we assume that the fast 
mode has been appropriately sampled and does not need to be down sampled 
since the user has control over the sampling rate. However, this result is eas- 
ily generalisable allowing both signals to be analysed at appropriately down 
sampled rates. 

In many non-linear systems there arc large numerical differences of the 
regressors due to the basis function(s) used to estimate the static map. The 
large differences lead to the regressor matrix being ill-conditioned and results 
in unstable matrix inversion and poor parameter estimates [11]. To alleviate 
ill-conditioning we propose using standardised least-squares (SLS) in combi- 
nation with the dual sampling approach proposed in Eqn. 5. 

Given a matrix of independent variables 0 and of dependent variables Z 
compute the mean and standard deviation of each variable, and replace 0 and 
Z with the centred and standardised variate as 


0 = (t-H*)- 


Z = (Z-p z )E z 1 


where E^ is a diagonal matrix of standard deviations with E^ denoting the 
standard deviation of the kth column of 0 and is a matrix who’s kth column 
has all entries equal to the mean of column k of 0. Substituting Eqn. 6 into 5 
yields a dual sampling rate SLS formulation in a sparse structure. 
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In the sequel we use the formulation in Eqn. 7 to estimate continuous-time 
parameters of nystagmus and saccade dynamics. 


4 Simulations & Results 

The accuracy of the dual sampling rate SLS technique with sparsity was vali- 
dated by simulating a continuous-time VOR model (e.g. Fig. 1) using Simulink. 
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Table 2. Left: Continuous-time coefficient values. T \ : slow-phase time- 
constant, t 2 : fast-phase time-constant, K \ : slow-phase gain, Kr. fast-phase 
gain and T: sampling interval. Right: Coefficient values of static non-linearity. 
b : lineal - term, c: squared term and d cubic term. 

The parameters used in the simulation are shown in Table 2 and representative 
of typical values found in experiments [4]. 

One hundred Monte-Carlo simulations were generated in which the input- 
output realisation was the same but had a unique Gaussian white, zero-mean, 
noise sequence added to the output. Excluding the noise free (NF) case the 
signal-to-noise ratio (SNR) of the noise sequence was decreased from 20 — 0 
dB in increments of 5 dB. The system was perturbed using a sinusoid input (1/6 
Hz frequency and 188 deg/s amplitude). A sinusoid input was used because it 
is the type of perturbation used in clinical settings. The system input-output 
was sampled at 100 Hz. 

The system parameters were estimated as outlined in Eqns. 4-7. Specifi- 
cally, for each input-output realisation, we analysed the VOR model as follows: 

1 . Single Sampling Rate, N on-Standardised Least-Squares : Both VOR modes 
were analysed using a single sampling rate and identified simultaneously 
using a sparse structure using traditional least-squares (Eqn. 4), 

2. Dual Sampling Rate, N on-Standardised Least-Squares : fast-phase was 
sampled at the original rate (100 Hz) whilst slow-phase was down sam- 
pled by 10(10 Hz) and identified simultaneously using a sparse structure 
using traditional least-squares (Eqn. 5), 

3. Dual Sampling Rate, Standardised Least-Squares: fast -phase was sam- 

pled at the original rate (100 Hz) whilst slow-phase was down sampled 
by 10 (10 Hz) and identified simultaneously using a sparse structure us- 
ing SLS (Eqn. 7). 

The continuous-time parameters of slow and fast-phase dynamics were esti- 
mated using the theoretical relationships in Table 1. Notice that we consider 
the lineal' system to have unity gain and translate the overall gain onto the poly- 
nomial basis function giving an estimate as a product of the linear system and 
static non-linearity gain. This is necessary because it is impossible to measure 
the signal at the output of the static non-linearity. Therefore, we deem that the 
best estimate of the linear system gain is a product of the linear system gain 
and lineal' coefficient of the static non-linearity, i.e. G\gb. 
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The results of this study arc shown in Fig. 2. The panels of the left col- 
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Figure 2. Left column: Single Sampling Rate, Non-Standardised Least- 
Squares. Centre column: Dual Sampling Rate, Non-Standardised Least- 
Squares (100 & 10 Hz). Right Column: Dual Sampling Rate, Standardised 
Least-Squares (100 & 10 Hz). Ordinate: STD about mean. Abscissa: Output 
SNR = NF, 20, 15, 10, 5, 0 dB, where NF denotes noise free case. (Note that 
the abscissa is shown in decreasing SNR which corresponds to increasing noise 
amplitude.) 

umn illustrate our findings for a non-standardised (traditional) least-squares 
approach using a single sampling rate to analyse both modes. The results show 
the nystagmus time constant is significantly biased with large variance whilst 
the nystagmus gain and saccade parameters are as theoretically expected with 
increasing bias and variance for decreasing SNR. The centre column displays 
the results of using a dual sampling rate with a non-standardised least-squares 
technique. The panels show this method improves bias and variance estimates 
for decreasing SNR since it yields theoretically expected trends. However, the 
non-monotonic increase in variance for the slow-phase time constant (e.g. 10-0 
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dB) suggests that the regressor matrix tends to becomes ill-conditioned for de- 
creasing SNR. The right column displays the results of using a dual sampling 
rate with a standardised least-squares approach. These results demonstrate that 
a dual sampling rate approach in combination with SLS significantly improves 
bias and variance estimates to expected trends. Nevertheless, the bias and vari- 
ance of the slow phase gain increase with dual sampling rate and standardis- 
ation. This result is expected since in these two cases the gain parameter is 
estimated with 10 times less data. Hence, using SLS the slow-phase variance 
monotonically increases for decreasing SNR and agrees with theoretically ex- 
pected results. 


5 Discussion 

Simulation results presented in §4 demonstrate that a dual sampling rate ap- 
proach provides improved results to single sampling rate techniques. In addi- 
tion, when combined with preconditioning offered by SLS, estimates are fur- 
ther improved. Hence, our dual sampling rate SLS technique is a robust and 
easily applicable methodology for the analysis of VOR dynamics. 

In many continuous and discrete-time parameter estimation problems the 
estimates may have incorrect sign due to high variance, possibly due to numer- 
ical ill-conditioning. Our proposed technique alleviates this problem by using 
matrix preconditioning (compare Fig. 2 (a)-(b) with (c)). This result may have 
general applicability for many estimation problems in biology such as for the 
analysis of ankle dynamics [12]. 

The central problem with VOR analysis is due to its time-constants (or 
poles) being separated by many orders of magnitude. In our example with a 
clinically relevant parameter set, the system poles arc separated by 30 orders of 
magnitude, which yields in slow-phase being highly oversampled and leading 
to numerical instability. 

The approach presented here is generalisable to most bifurcating systems 
since they have poles that arc typically several orders of magnitude apart, re- 
sulting in highly oversampled dynamics and unstable parameter estimates. 


6 Conclusions 

The results demonstrate that our dual sampling rate SLS approach with sparsity 
is a robust and easily implementable technique for the analysis or VOR data. 
These results suggest that clinical analysis of nystagmus (and saccade) may be 
improved, providing clinicians more accurate information for diagnosis, dis- 
ease modelling, etc. In addition, the robust estimation of VOR may lead to 
improved mission planning for astronauts and test-pilots of advanced aircraft. 
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